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Transport equation of the galactic cosmic ray (GCR) is numerically solved for qA > and qA < based 
on the stochastic differential equation (SDE) method. We have developed a fully time-dependent and three- 
dimensional code adapted for the wavy heliospheric current sheet (HCS). Results anticipated by the drift pattern 
are obtained for sample trajectories and distributions of arrival points at the heliospheric boundary for GCR 
protons. Our simulation reproduced a 22-year cycle of solar modulation which is qualitatively consistent with 
observations. Energy spectra of protons at 1 AU are calculated and compared with the observation by BESS. 



1. Introduction 

Transport of the galactic cosmic rays in the heliosphere is described by Parker's transport equation. The 
Parker's equation is equivalent to a coupled stochastic differential equation (SDE) [1, 2]. This SDE method 
allows us to get some information about solar modulation phenomena of galactic cosmic ray (GCR) not ob- 
tained by other numerical methods, such as distributions of arrival time, energy lost, and trajectory. We have 
developed a fully time-dependent and three-dimensional code based on the SDE method. We present some 
numerical results on solar modulation effects on GCR obtained by this new code. 



2. The Model 

The SDE equivalent to Parker's transport equation is written as (2l 



dP = \P{V-Y)dt, 



(1) 



where X and P are the position and the momentum of the particle, V is the solar wind velocity, Vd is the 
gradient-curvature drift velocity, n is the diffusion coefficient tensor, cr^cTg = 2k^'', and dWs is a Wiener 
process given by the Gaussian distribution. We adopted y = 400km/s, K|| — 1.5xlO^^/3(p/(lGey/c))(i?e/_B) 
cm^/s, Hi±_ — 0.05k||, B is the magnetic field, and Be ~ 5nT is B at the earth. We assume the structure of 
hehospheric magnetic field (HMF) as standard Parker spiral. We have used the "drift velocity field method" 
for the calculation of drift in the heliospheric current sheet (HCS). In our simulation, particles start at 1 AU 
on the equatorial plane and run backward in time until they exit the heliospheric boundary, 80 AU. 

We assume the structure of the HCS is as 
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where r, 6, (p are spherical polar coordinates relative to the solar rotation axis, (f>o is a constant, is the radius 
of the Sun, and Hq is the rotational angular velocity of the Sun. Similar expression for HCS had been adopted 
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Figure 1. (a) A view of the 
HCS projected onto a merid- 
ional plane as a function of 
radial distance from the Sun. 
The tilt angle near the Sun 
is 10°, and the solar activ- 
ity heading for the maximum; 
(b) Same as (a), but the solar 
activity heading for the mini- 
mum;(c) Same as (b), but the 
tilt angle near the Sun is 70°. 



by Jokipii and Thomas |4|, but with constant tilt angles contrary to our model. In this simple HCS model, 
the tilt angle changes continuously from 90° to 0° to 90° over a period of 11 year. The solar polarity of the 
HMF reverses when the tilt angle of HCS near the sun reaches 90°. Figure [2 shows views of HCS projected 
onto a meridional plane for three phases in the solar cycle. We have not taken into account in our model such 
phenomena as the global merged interaction regions, and the dynamically changing HCS with Fisk-like field 
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3. Results 

Figure |2] shows a sample trajectory for two protons which propa- 
gate from a point at 1 AU on the equatorial plane to 10 AU. The 
coordinate system of the figure is corotating with the sun. The 
drift pattern of the GCR depends on the polarity qA, where g is a 
charge sign of the sample particle and A is a constant character- 
izing HMF. For the positive polarity qA > 0, the particles drift 
from the polar region towards the HCS, and outward along the 
HCS. In contrast, for the negative polarity qA < 0, the direction 
of the drift is opposite to the case for qA > 0. 

Figure |3] shows the heliolatitude-longitude distribution of simu- 
lated particles when they exit the heliospheric boundary, where 
the longitude is shown in the corotating coordinate system. The 
thick line indicates the HCS at 1 AU when the particles start off 
from a point (indicated by "star") at 1 AU on the equatorial plane. 
Latitudes of the arrival points for qA < distribute around the 

HCS as expected, because particles arrive at 1 AU along the HCS. The thickness of the distribution of latitude 
(±30°) reflects the structure of spatiotemporal variation of the HCS. The structure of the HCS seen by particles 
is different from particle to particle depending on the arrival time. On the contrary, the latitude for qA > dis- 
tribute around the polar region, either northern (a) or southern (b) depending on whether the observer is located 
at a point northern (a) the HCS or southern (b), respectively. This tendency is also as expected, because it is 
difficult for positive particles to penetrate through the HCS during the journey from the outer heliosphere. The 
mean arrival time of 1 GeV protons is 4.5 days and 80 days for qA > and qA < 0, respectively. The short 



Figure 2. Sample trajectory for two protons. 
The tilt angle near the Sun is 40°, and the solar 
activity heading for the maximum. The energy 
at 1 AU on the equatorial plane is 1 GeV. 
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Figure 3. (a) Heliolatitude-longitude distribution by simulation for protons witli energy 1 GeV at 1 AU when they exit the 
heliospheric boundary, 80 AU. The observer at 1 AU on the equatorial plane is located at a point northern the HCS; (b) 
Same as (a), but the observer is located at a point southern the HCS. The tilt angle near the sun is 40°, and the solar activity 
heading for the maximum. 



time of only 4.5 days for qA > suggests that particles are less modulated. This relatively little modulation 
effect is reflected in the fail to reproduce the observed spectrum for qA > as we mention later. 

Figure]^ shows simulated variation of the proton intensities over the 22-year modulation cycle. The solid line 
indicates the result for the model of the HCS by eq|2] («(?', t) model). The dashed line shows the result for 
the HCS model (a{t) model) where the tilt angle varies with time but does not depend on radius from the sun. 
In our a{r, t) model, the variation of the HMF near the sun propagates with the solar wind velocity. It takes 
about 1 year for the signal of polarity change and change in the tilt angle reach the heliospheric boundary. In 
a{t) model, on the contrary, these changes reach instantaneously the boundary. This difference between the 
two models in signal propagation is reflected in a time lag of the proton intensity changes by a{r, t) model as 
seen in Figure0t. Our model succeeded in reproducing a 22-year modulation cycle and in reproducing, though 
qualitatively, rather flat maximum for qA > and a sharp maximum for qA < . 

Figure IIJ) shows proton energy spectra by our simulation (solid line) and the spectra observed by BESS ex- 
periment (open symbol). The open triangles refer to proton fluxes observed in 1998 |7| when the tilt angle 
near the sun is about 10° and the solar polarity is positive. The open squares refer to proton fluxes observed 
in 2000 1 8 1 when the tilt angle near the sun is about 70° and the solar polarity is negative. The symbols for 
simulated energy spectra refer to the phases marked by the same symbol in Figure]^. The view of the HCS 
corresponding to these phases are shown in Figure^ andfO;. The local interstellar spectrum (LIS) of proton is 
assumed as jr oc /3(r + O.Smc^)^^'^, where T is kinetic energy and m is the mass of proton. The simulated 
energy spectra are normalized at 3.6 GeV to the observation in 2000 |8|. Our model spectrum for qA < 
agrees quite well with BESS result, however our model fails to reproduce the observed spectrum for qA > 0. 
The failure of our model for qA > may come from inadequacies of the heliospheric structure in our model. 
As we mentioned before, positive particles drift from polar region within short times modulated rather weakly. 
If we take into account the existence of the random transverse component of the HMF in the polar regions ||9l, 
we may be able to reproduce the observed spectrum because it is expected that the random transverse compo- 
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Figure 4. (a) Simulated variation of ttie 1 GeV proton intensities over the 22-year modulation cycle; (b) Simulated proton 
energy spectra and observed proton energy spectra in 1998 and 2000. 

nent suppresses drifts in the polar regions. The actual variation of the HCS is much irregular in contrast to our 
model. We have to consider the actual dynamically changing HCS. 

4. Conclusions 

We have developed a time-dependent and three-dimensional code based on the SDE method. The sample 
trajectories and distributions of arrival points at the heliospheric boundary agree with the physical expectation. 
Our simulation reproduced a 22-year modulation cycle which is qualitatively consistent with observations. 
Our model spectrum for qA < Q agrees with BESS result, however our model fails to reproduce the observed 
spectrum for qA> 0. 
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